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Recent  advances  in  computational  techniques  have  allowed  the  application  of  computational  tools  to 
study  heterogeneous  functional  materials  (HeteroFoaMs)  in  the  solid  oxide  fuel  cell  (SOFC)  from  the 
quantum  (sub-atomic)  to  atomistic  to  the  continuum  scales.  However,  knowledge  gained  from  a  par¬ 
ticular  computational  technique  can  only  provide  insight  at  that  specific  scale.  There  has  been  a  recent 
interest  to  develop  a  more  cohesive  effort  so  that  results  obtained  from  models  across  a  particular  spatial 
dimension  can  be  used  to  extract  additional  insight  across  a  larger  range  of  length  scales.  This  review 
article  surveys  recent  progress  in  the  modeling  and  simulation  of  SOFCs,  and  relates  them  to  the  rel¬ 
evant  physical  phenomena  and  length/time  scales.  We  then  proceed  to  review  the  various  numerical 
techniques  used,  and  their  applicability  across  the  length  and  time  scales. 

©  201 1  Elsevier  B.V.  All  rights  reserved. 
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1.  Background  and  introduction 

Modeling  and  simulation  provides  a  unique  opportunity  to 
assist  with  the  development  of  materials,  components,  and  prac¬ 
tices  for  next  generation  solid  oxide  fuel  cells  (SOFCs).  The  methods 
enhance  our  ability  to  comprehend  the  connections  between  func¬ 
tional  and  detrimental  behavior  and  their  impact  on  the  cell. 
However,  it  is  a  challenge  to  treat  these  aspects  across  the  breadth 
of  time  and  length  scales  over  which  they  originate  and  manifest 
themselves.  This  is  complicated  by  the  heterogeneous  nature  of  the 
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material  structures,  which  convolutes  the  challenge  of  understand¬ 
ing  these  processes.  Modeling  and  simulation  tools  can  help  us 
to  understand  these  complex  interactions.  With  continued  devel¬ 
opment,  these  methods  may  eventually  be  able  to  be  applied  to 
provide  a  predictive  capability  for  the  design  and  development 
of  SOFC  materials,  morphologies/structures,  components  and  sys¬ 
tems;  however,  this  remains  a  significant  challenge  [1  ]. 

As  with  the  heterogeneous  functional  materials  (HeteroFoaMs) 
[167]  themselves,  the  individual  constituent  materials,  interfaces, 
and  morphologies  influence  the  approach(es)  taken.  The  physical 
and  (electro)chemical  conditions  of  the  materials  and  processes 
involved  can  limit  the  validity  and  range  of  applicable  approaches. 
The  relevant  time  and  length  scales  compound  the  challenge  as  they 
ultimately  dictate  the  plausible  approaches  that  may  be  taken.  In 
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Fig.  1.  An  overview  of  the  time  and  length  scales  characteristic  of  functional  and  degradation  processes  in  the  heterogeneous  SOFC  electrode  structures.  The  time  and  length 
scales  provided  are  conceptual  and  are  only  intended  to  be  characteristic  of  the  scales  at  which  these  processes  generally  present  themselves. 


Fig.  1,  time  and  length  scales  representative  of  SOFC-relevant  func¬ 
tional  and  degradation  processes  are  provided.  The  interactions  for 
each  of  these  phenomenon  originates  at  the  quantum  level  with 
the  electronic  structure  of  the  individual  atoms;  however,  they 
present  themselves  in  unique  observable/measurable  ways  which 
can  occur  at  a  variety  of  time  and  length  scales  and  also  provide 
insights. 

Examining  Fig.  1,  chemical  and  electrochemical  processes  are 
governed  by  fundamental  aspects  of  the  catalytic  and  electrocat- 
alytic  interface.  These  processes  are  dictated  by  the  chemical  and 
physical  interactions  of  the  molecules  and  catalytic  binding  sites 
to  Refs.  [2-6].  These  processes  present  themselves  at  the  time 
and  length  scales  of  the  electronic  interactions  of  the  individ¬ 
ual  atoms/molecules  that  determine  bond  strengths  and  binding 
energies  of  the  adsorbed  species  on  up  to  those  of  the  transport 
processes  which  supply  them.  Fracture  and  stresses/strains  involve 
the  interactions  at  the  lattice  level  and  are  directly  influenced  by 
features  including  grain  interfaces,  phase  interfaces  and  defects  on 
up  through  the  microstructure  [7-11]. 

Homogeneous  transport  includes  mass,  charge,  and  heat  trans¬ 
fer,  which  occurs  through  the  bulk  of  the  pores,  materials,  and 
structures.  These  processes  ranges  from  time  and  length  scales 
of  the  materials  lattice  and  structures/features  that  support  them 
(at  the  corresponding  thermodynamic  state)  on  up  to  the  systems 
levels  [12-15].  On  the  other  hand,  surface  and  grain  bound¬ 
ary  transport  processes  must  match  the  length  and  time  scales 
attributed  to  the  grain  structure  itself  [16-19]. 

Aging  (e.g.,  coarsening),  sintering,  and  redox  cycling  tend  to 
occur  at  scales  attributed  to  the  grain  structure  and  microstructure 
as  they  involve  chemical  and  physical  changes  in  the  materials, 
interface,  and  morphology  [20-26].  Similarly,  poisoning,  coking, 
and  passivation  involve  changes  to  the  materials  and  interfaces 
through  chemical  interactions  meaning  these  aspects  also  present 
themselves  at  the  scales  of  the  lattice  and  grain  structure  and  can 


impact  broader  scales  [20,27-29].  At  the  broader  time  and  length 
scales  are  aspects  like  the  bulk  mechanical  properties  and  compo¬ 
nent/systems  responses  to  the  conditions  and  history  (e.g.,  creep 
and  fatigue)  [11,30-33].  These  processes  include  more  fundamen¬ 
tal  aspects,  such  as  the  coefficients  of  thermal  expansion,  but  are 
directed  toward  the  bulk  components/cells. 

The  modeling  and  simulation  methods  that  are  applied  require 
comparable  levels  of  precision  and  control,  while  also  being  able  to 
accurately  treat  the  relevant  processes  and  conditions.  For  exam¬ 
ple,  specifics  of  the  mechanisms  attributed  to  the  electrochemical 
reaction  kinetics  may  be  studied  by  using  density  functional  the¬ 
ory  (DFT).  Various  processes  and  pathway(s)  can  be  determined  by 
mapping  the  corresponding  potential  energy  surfaces,  with  saddle 
points  corresponding  to  transition  states.  However,  semiconduc¬ 
tors  and  materials  with  highly  correlated  electronic  structures 
require  more  rigorous  approaches,  such  as  DFT  +  U  and  quantum 
Monte  Carlo  methods,  because  the  exchange  correlation  functions 
of  traditional  DFT  methods  do  not  provide  the  necessary  accu¬ 
racy  [1,34].  This  can  be  relevant  to  transition  metal  oxides  and 
ceramics,  like  those  often  used  in  the  SOFC.  Aspects  requiring 
the  incorporation  of  longer  range  and  non-periodic  effects  can 
also  be  difficult,  especially  as  finite  size  effects  including  hetero¬ 
geneities  and  morphologies  come  into  play.  Processes  involving 
chemisorption,  vacancies/defects  or  impurities,  and  electrochemi¬ 
cal  polarization  are  examples,  which  generally  require  embedded 
cluster  methods.  Accurate  partitions  are  required  to  couple  phe¬ 
nomena,  which  are  difficult  and  specific  problems  in  their  own 
right;  however,  again  necessary  [34].  Insufficient  treatment  of  these 
aspects  can  result  in  un-reliable  predictions  in  mechanisms  and 
properties  attributed  to  the  system. 

At  broader  scales,  there  are  parallels  in  the  methods  for  the  mod¬ 
eling  and  simulation  of  transport  phenomena,  where  aspects  such 
as  the  inclusion  of  the  appropriate  physical  phenomena  when  the 
length  scales  approach  the  mean  free  path  [13-15],  space  charge 
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Fig.  2.  Several  of  the  modeling  and  simulation  methods  applicable  for  the  SOFC,  where  the  scalable  cost  is  plotted  versus  the  time  and  length  scales  that  can  be  resolved. 
Computational  cost  represents  relative  trends. 


effects  [35-38],  the  presence  of  poisons  [27,39-54],  and  so  on  can 
have  similar  influences.  Scenarios  like  the  presence  of  unexpected 
secondary  trace  elements  and  poisons  signify  a  limitation;  theoret¬ 
ical  approaches  require  an  accurate  understanding  of  the  physical 
conditions  germane  to  the  problem  at  hand.  If  unknown  traces 
of  poisoning  elements/species  leach  into  the  cell  from  the  sur¬ 
rounding  environment  and/or  balance  of  plant,  their  contributions 
and  effects  cannot  be  accounted  for.  This  highlights  the  continued 
importance  and  need  for  detailed  experimental  efforts  for  checking 
and  validating  our  understanding  as  well  as  defining  the  constraints 
of  our  physical  problem. 

2.  Approaches  to  modeling  and  simulation 

The  typical  scientific  approach  to  research  and  development 
either  relies  upon  trial  and  error,  the  so  called  Edisonian  approach, 
or  it  relies  upon  the  maturation  of  the  scientific  understanding 
of  the  relevant  materials  properties,  thermodynamic  properties, 
physics,  chemistries,  and  responses  of  the  system  through  system¬ 
atic  exploration.  Individual  aspects  of  the  system  are  separated 
and  studied  to  help  refine  this  understanding.  This  has  often  been 
accomplished  by  means  of  tedious  experimental  effort.  Model¬ 
ing  and  simulation  across  the  time  and  length  scales  provide 
unique  opportunities  to  understand  properties  and  phenomena 
that  are  inaccessible  and  difficult  to  address  with  experimental 
effort,  or  have  inseparable  multifaceted  interactions.  As  noted  in 
Fig.  1,  discrepancies  in  the  length  and  time  scales  of  these  func¬ 
tional  and  degradation  processes  exemplify  these  modeling  and 
simulation  needs  for  understanding  the  porous  structure  and  its 
constituent  materials,  interfaces,  and  morphology.  Many  of  the  typ¬ 
ical  approaches  are  described  in  Fig.  2,  which  presents  the  modeling 
or  simulation  methods’  relative  computational  cost  (not  to  scale) 
versus  the  approximate  characteristic  time  and  length  scales  that 
they  can  resolve.  Comparing  Fig.  1  to  Fig.  2,  it  is  apparent  that  no 
single  method  is  suited  to  studying  the  full  scope  of  challenges. 


Fig.  2  demonstrates  that  different  approaches  are  required  to 
capture  different  time  and  length  scales.  While  aspects  including 
the  relative  domain  size  and  durations  can  affect  the  computational 
cost  of  a  study,  the  cost  tends  to  decrease  with  increasing  time  and 
length  scales.  This  is  a  general  observation  that  accounts  for  the 
scalability  of  the  approaches  (i.e.,  parallelization),  but  should  not 
be  taken  as  a  universal  rule.  At  the  most  fundamental  scales  are 
ab  initio  Molecular  Dynamics  (AIMD)  and  ab  initio  methods,  such 
as  DFT,  examine  the  electronic  structure  (quantum)  of  the  respec¬ 
tive  atoms  and  the  corresponding  interactions  [1,34,55,56].  AIMD 
uses  the  resulting  potentials  and  forces  to  calculate  the  dynam¬ 
ics  of  the  atomic  and  molecular  interactions.  The  expense  of  these 
methods  is  tied  to  the  time  scales  of  the  electronic  structure  of  the 
atoms  relative  to  those  of  the  atomic/molecular  dynamics,  as  well 
as  the  numerous  interactions  between  atoms  and  molecules.  To  a 
lesser  degree,  this  holds  for  Molecular  Dynamics  and  Monte  Carlo 
approaches,  which  use  predetermined  force  fields  to  dictate  molec¬ 
ular  interactions  in  place  of  detailed  calculations  of  the  electronic 
structure  [56]. 

Discrete  element  and  phase  field  methods  are  often  used  to 
study  the  formation  and  evolution  of  grain  and  microstructure 
in  such  systems.  This  may  include  aspects  such  as  sintering, 
compaction,  and  aging  resulting  from  the  non-equilibrium  ther¬ 
modynamics  in  such  systems  [57-59].  Transport  phenomena  in 
discrete  microstructures,  may  be  captured  using  methods  built 
upon  the  Boltzmann  transport  equation,  such  as  the  lattice  Boltz¬ 
mann  method  (LBM),  which  is  scalable  and  enables  the  use  of 
a  low  resolution,  regular  mesh  [60-62].  The  underlying  Boltz¬ 
mann  equation  is  applicable  from  free  molecular  to  the  continuum 
scales;  however,  LBM  is  typically  formulated  so  that  it  recovers  the 
Navier-Stokes  and  energy  equation.  Empirical  and  rigorous  mathe¬ 
matical  approaches  have  been  attempted  to  extend  its  applicability 
to  finite  Knudsen  numbers  [63-65].  More  traditional  methods, 
such  as  those  using  the  Navier-Stokes  equations,  energy  equation, 
dusty  gas  equations,  as  well  as  Navier’s  and  related  mechanical 
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constitutive  equations  may  also  be  considered  [  1 2-1 5,66,67].  These 
methods  are  typically  discretized  with  finite  element  and  finite 
volume  method,  which  can  be  difficult  to  mesh  in  complex  struc¬ 
tures.  Therefore,  these  methods  are  often  used  for  volume  averaged 
representations  of  the  materials/system  (i.e.,  homogenized).  The 
typical  numeric  methods  used  to  solve  these  equations  (i.e.,  finite 
element,  finite  difference,  and  finite  volume)  may  also  be  used  to 
solve  other  interactions  and  processes  as  long  as  they  are  prop¬ 
erly  represented  in  the  form  of  systems  of  partial  differential 
equations. 

Chemical  and  electrochemical  kinetics  are  typically  implicitly 
included  in  atomistic  and  molecular  methods;  however,  they  must 
be  determined  and  explicitly  accounted  for  at  broader  time  and 
length  scales.  Both  heterogeneous  and  homogeneous  kinetics  are 
typically  described  as  being  stiff,  indicating  a  disparity  of  time 
scales  involved  [4-6].  While  the  underlying  interactions  follow 
quantum  to  molecular  time  and  length  scales,  the  equilibration 
process  can  follow  and/or  dictate  transport  processes  indicat¬ 
ing  that  they  can  change  over  a  wide  range  of  time  and  length 
scales.  Solution  of  such  a  mechanism  is  repetitive,  can  dictate  time 
steps  permissible  for  the  modeling  of  other  coupled  processes,  and 
become  prohibitively  expensive  to  evaluate  as  the  complexity  of 
the  mechanism  increase. 


2.1.  System ,  cell  and  electrode  level  approaches 

Systems,  single-cell,  and  electrode-level  SOFC  modeling  con¬ 
stitute  one  of  the  broadest  efforts  to  represent  the  SOFC  and  its 
function  with  theory.  These  approaches  are  often  broken  down 
into  steady  state  or  transient  models  and  can  range  from  zero 
to  three  dimensional  descriptions  of  the  system,  cell  and/or  elec¬ 
trode.  These  approaches  are  often  used  to  study  the  effects  of 
changes  in  component  design/materials,  fuel  and  oxidant,  condi¬ 
tions,  and/or  coupling  between  processes.  Specific  aspects  of  the 
system,  cell  or  electrode,  such  as  cell  performance,  fuel  utilization, 
stresses/strains,  dynamic  response,  and/or  thermal  signatures  and 
conditions  are  typically  the  focus  of  such  efforts. 

To  describe  interactions  in  such  a  system,  conservative  contin¬ 
uum  theory,  such  as  Ohm’s  law,  Fields  law  or  a  Stefan-Maxwell 
or  Onsager’s  based  multi-component  theory,  Navier-Stokes  equa¬ 
tions,  energy  equation,  Navier’s  equation  or  a  related  small 
deformation  continuum  mechanics  theory,  is  often  used  with 
the  inclusion  of  non-continuum  principles  as  applicable.  Non¬ 
continuum  principles  can  become  important  for  aspects  including 
the  description  of  diffusive  mass  transfer  because  of  the  tempera¬ 
ture  and  length  scales  involved  (e.g.,  Knudsen  flow).  These  aspects 
are  often  treated  with  the  use  of  the  dusty  gas  equations  [15].  Cell 
potentials  are  often  dictated  by  the  Nernst  equation  with  activa¬ 
tion  overpotentials  calculated  via  the  Butler-Volmer  equation.  In 
general,  finite  elements,  finite  differences,  finite  volumes,  colloca¬ 
tion  methods,  and/or  marching  schemes  are  used  to  discretize  and 
solve  the  governing  equations. 

This  is  representative  of  the  three-dimensional  models  pub¬ 
lished  by  Khaleel  and  coworkers,  which  couple  transport  and 
electrochemical  processes  within  SOFC  single-cells  and  stacks 
[68-70].  Finite  element  analysis  was  used  to  solve  a  Navier-Stoke’s 
based  description  flow,  along  with  the  energy  equation  to  account 
for  heat  transfer  in  the  form  of  conduction  and  advection  [68,69]. 
An  electrochemical  model  with  the  Butler-Volmer  equations  was 
used  to  calculate  the  overpotentials  associated  with  the  oxida¬ 
tion  and  reduction  reactions,  along  with  the  corresponding  heats 
of  reaction.  The  effects  of  fuel  and  oxidant  flow  configurations 
and  cell  layout  on  dynamic  and  localized  temperature  distribu¬ 
tions,  concentrations,  and  pressures  have  been  examined  with  this 
approach. 


A  key  consideration  for  this  approach  is  how  the  het¬ 
erogeneous  electrode  structure  should  be  represented  and  its 
influence  on  the  functional  processes.  Empirical  factors,  such  as 
a  porosity-tortuosity  factor,  exchange  current  density,  TPB  length, 
catalytic  area,  percolation  probability  or  effective  media  approx¬ 
imations,  permeability  coefficients,  mean  pore/particle  size,  and 
area  specific  resistances,  are  used.  These  parameters  may  be 
obtained  using  a  variety  of  approaches  including  fitting  and  esti¬ 
mating  from  independent  experimental  measurements.  As  will  be 
discussed  in  more  detail,  theoretical  and  conceptual  representa¬ 
tions  of  the  microstructure  are  also  commonly  used  to  estimate 
these  types  of  parameters.  The  underlying  premise  of  this  approach 
is  to  simplify  the  conceptual  construct  so  that  the  cell’s  physical 
behavior  and  responses  may  be  replicated.  Then,  a  full  gamut  of 
variations  occurring  within  the  cell  can  then  be  obtained  and  the 
model  may  be  validated  using  experimental  data.  Parametric  vari¬ 
ations  in  the  physical  cell  description  or  the  empirical  parameters 
can  then  be  used  to  optimize  aspect(s)  of  the  cell  design. 

The  importance  of  the  microstructure’s  description  has  been 
noted  by  Kee  and  coworkers  [71-73].  Cell-level  modeling  efforts 
provide  a  valuable  tool  for  obtaining  these  types  of  details;  how¬ 
ever,  making  adjustments  to  the  various  cell  components  based  on 
these  global  models  may  be  challenging.  The  methods  implicitly 
assume  that  the  complete  sets  of  appropriate  physics  are  being  con¬ 
sidered  with  sufficient  fidelity  to  capture  the  spectrum  of  relevant 
processes.  Complex  heterogeneous  structures  are  treated  on  a  vol¬ 
ume  averaged  (homogenized)  basis,  assuming  consistency.  While 
this  provides  an  effective  means  of  making  it  tractable  to  solve, 
disconnects  can  exists  between  the  identification  of  target  values 
for  these  empirical  parameters  in  these  models  and  how  the  het¬ 
erogeneous  structure  may  be  changed  to  achieve  the  target.  More 
specifically,  how  this  may  be  done  without  having  other  deleterious 
effects. 

An  interesting  and  unique  attempt  to  address  some  of  these 
challenges  has  been  taken  up  by  Bessler  and  coworkers  [74-77]. 
Transport  and  reaction  models  that  calculate  electrochemical 
impedance  spectra  have  been  developed  so  that  electrochemical 
and  transport  contributions  incorporated  in  the  cell’s  description 
can  be  translated  into  the  frequency  domain.  This  allows  con¬ 
tributions  within  the  physical  and  electrochemical  processes,  as 
well  as  those  related  to  the  description  of  the  microstructure, 
to  be  separated  based  upon  their  contribution  in  the  frequency 
domain.  Comparing  to  complementary  experimental  electrochem¬ 
ical  impedance  spectroscopy  data  provides  confidence  in  the  details 
extracted.  The  details  ascertained  from  the  impedance  model¬ 
ing  can  then  be  used  to  simulate  cell  transport  and  polarizations 
[77,78]. 

A  rather  comprehensive  cell-level  approach  has  been  reported 
by  Kee  and  coworkers,  which  involves  the  direct  internal  hetero¬ 
geneous  methane  reformation  reaction  kinetics  developed  and 
incorporated  into  two-  and  three-dimensional  transport  models 
of  planar  and  tubular  SOFCs  [71,79-82].  Elementary  reaction 
kinetics  describing  internal  methane  reformation  processes  have 
been  assembled  for  this  effort  and  are  numerically  solved  with  a 
damped-Newton  method  for  nonlinear  systems  of  equations  using 
the  pseudo-steady  state  approximation.  Mass  transport  in  the 
porous  electrode  is  treated  with  the  dusty  gas  equations,  which 
is  coupled  with  principles  of  species,  electronic  and  ionic  charge, 
continuity,  momentum,  energy  conservation  within  the  system 
[71,80-86].  An  extensive  effort  has  been  placed  on  validation  of 
their  reaction  kinetics  model  using  custom  flow  cell  with  a  isolated 
Ni-YSZ  anode  and  mass  spectroscopy  probes  to  measure  changes 
in  chemical  compositions,  which  showed  their  kinetics  were  able 
to  reproduce  the  experimental  measurements  with  considerable 
fidelity  over  a  wide  range  of  flow,  composition  and  thermal 
conditions  corresponding  to  SOFC  operation  [87].  Additionally, 
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direct  internal  reformation  cell  models  exhibited  agreement  with 
SOFC  polarization  experiments,  including  for  various  methane  feed 
compositions  and  conditions  [71,80-86].  Results  from  these  works 
have  provided  insight  into  the  nature  and  extent  of  heterogeneous 
internal  methane  reformation  in  the  SOFC.  In  collaborative  efforts 
with  Barnett  et  al.,  this  approach  has  been  extended  to  study  coking 
processes  [86].  This  has  helped  to  spurn  design  developments, 
such  as  anode-oxide  support  layers  to  enhance  chemical  stability 
with  hydrocarbon-based  fuels,  where  the  models  were  able  to 
provide  the  proof  of  concept  and  design  input  [88]. 

Among  the  challenges  confronting  the  heterogeneous  SOFC 
structures  is  to  maintain  mechanical  integrity.  Elevated  temper¬ 
atures  and  corresponding  gradients  cause  the  materials  to  undergo 
thermal  expansion.  Differences  in  coefficients  of  thermal  expansion 
(CTE)  between  the  constituent  materials  and  layers  of  the  cell  result 
in  thermal  stresses.  Thermal  stresses  can  cause  cracking  within  the 
microstructure  or  even  delamination  and  other  catastrophic  fail¬ 
ure  [  1 1 ,30,32,89-91  ].  These  effects  are  most  pronounced  during  the 
transient  period  associated  with  thermal  cycling  and  load  following 
that  can  induce  large  internal  temperature  gradients  [31,32,92,93]. 
They  may  be  further  compounded  by  stresses  induced  by  the 
mechanical  constraints  from  the  cell’s  seals  and  packaging  [31,94], 
volumetric  changes  during  redox  cycling  [32,95]  and  poisoning 
[54]. 

In  many  regards,  efforts  to  simulate  the  thermal  stresses  at  the 
system,  cell,  and  electrode-levels  are  similar  and  complementary  to 
the  performance  models.  These  models  require  that  the  materials 
and  structures  be  treated  on  an  effective,  homogenized  basis.  Finite 
element  analysis  is  often  used  to  solve  the  thermo-elastic  prob¬ 
lem  (e.g.,  Energy  and  Navier’s  equations)  along  with  the  required 
constitutive  equations,  which  assume  a  linear  material  response 
according  to  the  principles  of  small  and  elastic  deformations 
[30,32,90,93,94,96].  Challenges  to  these  approaches  include  the  (i) 
treatment  of  the  often  brittle  material  as  elastic  or  quasi-elastic, 
(ii)  determination  of  the  baseline  mechanical  properties  which  can 
be  strongly  dependent  upon  microstructure  and  composition,  (iii) 
development,  coupling,  and  validation  of  a  detailed  thermal-fluid 
and  chemical  models,  (iv)  identification  and  implementation  of 
appropriate  boundary  conditions,  and  (v)  validation  of  the  models. 
Validation  can  be  difficult  because  system  and  operational  aspects 
play  a  prominent  role;  however,  methods  including  nanoindenta¬ 
tion  [33]  and  X-ray  diffraction  (XRD)  [97]  have  provided  some  new 
capabilities. 

Attempts  to  extend  the  stresses/strain  fields  calculated  by  the 
thermo-mechanical  and  degradation  models  have  also  been  used  to 
predict  crack  initiation  and  propagation  dynamics  [10,11,30,89,92]. 
Crack  propagation  modes,  trajectories,  stress  intensity  factors,  and 
interactions  resulting  from  residual  stress  and  strain  fields  must 
be  considered  [10,89].  The  intermittent  and  unpredictable  nature 
makes  it  a  formidable  challenge  to  assign  parameters  and  vali¬ 
date  results.  Weibull  statistics  have  been  incorporated  to  predict 
failure  probabilities  resulting  from  mechanical  degradation  mecha¬ 
nisms  [11,30-32,90,92].  A  probabilistic  measure  of  failure  provides 
a  measure  of  the  influence  of  different  aspects  of  the  heteroge¬ 
neous  structure  and  cell  design.  These  principles  have  been  used 
to  demonstrate  limits  on  thermal  gradients  [32],  cycle  steps  and 
conditions  [32],  contact  pressure  and  flow  conditions  [31,32],  and 
functional  grading  of  the  electrode  structure  [90].  Flowever  insight¬ 
ful,  these  analyses  can  be  difficult  to  validate  without  including 
far-reaching  experimental  programs. 

This  sample  perspective  of  the  system,  cell,  and  electrode  level 
models  that  have  been  reported  within  the  SOFC  community  pro¬ 
vides  detail  on  the  nature  and  impact  that  such  approaches  can 
have.  At  these  time  and  length  scales,  more  extensive  reviews 
on  the  dynamics  [98],  aspects  of  operation  [99],  and  modeling 
approaches  [99,100]  are  also  widely  available. 


2.2.  Discrete  morphological  models 

We  move  from  cell  level  to  discrete  model  representations  of 
the  heterogeneous  structure’s  morphology.  This  class  of  models 
are  described  as  morphological,  meaning  that  they  consider  either 
an  (i)  discrete  idealized  or  (ii)  artificially  generate  representation 
of  the  heterogeneous  structure.  These  systems  are  used  to  test 
microstructure  properties  and/or  to  test  hypothesis  regarding  the 
effects  of  properties  selected  during  manufacturing. 

Discrete  idealized  structures  may  be  considered  as  those  where 
physical  packing  of  the  microstructure  is  fixed.  Percolation  theory, 
which  uses  the  coordination  of  electrode  and  electrolyte  “parti¬ 
cles”  in  a  presumed  lattice  arrangement,  is  consistent  with  this 
definition  [73,101-104].  A  contact  angle  that  dictates  the  con¬ 
stricted  contact  between  discrete  particles  is  selected  as  an  input 
to  these  models.  However  idealized,  the  model  provides  an  esti¬ 
mate  of  the  reduction  in  conductivity  due  to  the  organization  of 
the  homogenized  microstructure.  The  effect  of  the  structure  on 
the  overall  conductivity  can  be  calculated,  along  with  estimations 
of  the  electrochemically  active  three-phase  boundary  (TPB)  length 
and  interfacial  areas  between  constitutive  particles  on  the  basis  of 
coordination  number,  particle  sizes,  volume  fractions,  and  average 
inter-particle  contact  angle. 

Results  from  percolation  theory  approaches  include  limitations 
on  volume  fractions  and  the  contributions  of  unique  particle  sizes 
to  TPB  length  and  degree  of  percolation  for  electronic  charge,  ionic 
charge  and  mass  transfer  processes  definition  [73,101-104].  Many 
of  the  studies  discussed  use  percolation  theory  to  assign  proper¬ 
ties  to  the  microstructure.  Percolation  theory  is  widely  accepted 
because  of  its  ability  to  reduce  the  complex  heterogeneous  con¬ 
struct  to  a  few  descriptive  parameters  that  follow  established 
trends,  including  validation  versus  experimental  area  specific  resis¬ 
tance  measurements  (ASR)  for  unique  electrode  compositions 
[101,104];  however,  it  is  inherently  limited  to  describing  these 
structures  as  highly  idealized  constructs.  Bimodal  particle  distribu¬ 
tions,  commonly  used  with  the  electrolytic  material  in  the  electrode 
structures,  are  difficult  to  treat.  Even  though  microstructure  prop¬ 
erties  can  be  estimated,  there  may  still  be  significant  variations 
within  actual  structures. 

A  second  idealized  approach  was  reported  by  Virkar  and 
coworkers  [105,106].  Rather  than  using  discrete  coordinated  par¬ 
ticles,  Virkar  and  coworkers  considered  electrolytic  fins  with 
electrode  particles  placed  on  the  periphery  of  the  fin  to  provide 
the  electrochemically  active  interfaces.  The  electrode  structures 
were  treated  as  homogenized  media,  with  effective  properties  that 
are  descriptive  of  the  microstructure.  These  fins  were  then  used 
to  calculate  effective  charge  transfer  resistances  on  the  basis  of  an 
ideal  intrinsic  charge  transfer  resistance  and  descriptions  of  the 
microstructure.  Effective  charge  transfer  resistances  were  found 
to  be  higher  for  thin  electrodes  and  then  decrease  once  a  criti¬ 
cal  electrode  length  was  reached.  Finer  particles  also  decreased 
effective  charge  transfer  resistances  [105,106].  In  these  studies,  the 
effective  charge  transfer  resistance  calculations  were  compared  to 
corresponding  experimental  measurements  where  the  electrode 
thickness  was  varied. 

In  a  similar  approach,  Nelson  et  al.  treated  the  structure  as 
an  electrochemical  fin  to  better  understand  charge  transport 
within  a  conceptual  structure  with  an  allowance  for  variable  cross 
sections  [107].  While  the  approach  includes  empirical  treatments 
of  particle  size  and  details  regarding  the  sinter  quality,  it  enabled 
distinct  charge  transport  regimes  within  the  SOFC  electrode  to  be 
understood  in  terms  of  how  they  relate  back  to  the  microstructure. 
The  study  was  validated  by  comparing  polarization  resistances 
for  both  electrodes  of  different  thicknesses  as  well  as  those  that 
were  well  and  poorly  sintered,  resulting  in  different  aspects  of 
inter-particle  contact.  An  interesting  result  of  this  approach  was 
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the  analog  to  traditional  heat  transfer  “fin  efficiencies”  which 
provide  an  effective  measure  of  the  quality  and  performance  of 
the  structure  relative  to  if  transport  and  electrochemical  reaction. 

An  alternative  approach  is  to  artificially  generate  a  structure  that 
may  be  used  to  represent  a  heterogeneous  material  that  may  be  rep¬ 
resentative  of  the  heterogeneous  electrode  structure.  These  types 
of  methods  rely  on  stochastic  procedures.  In  works  by  Sunde,  Monte 
Carlo  methods  were  used  [  1 08-1 10].  Two  approaches  were  shown 
by  Sunde.  The  first  used  random  number  generators  to  populate 
a  cubic  lattice  with  electrode,  electrolyte,  and  pore  particles  with 
some  pre-determined  volume  fraction  [108].  In  a  second  more  rig¬ 
orous  approach,  the  random  placement  of  spherical  particles  was 
used  [109,1 10].  Additional  particles  were  randomly  brought  in  and 
added  to  a  three  particle  cluster  of  spheres  in  hard  contact.  Particle 
contact  was  fixed  as  a  discrete  point.  Overlap  for  both  approaches, 
representative  of  the  sinter  process,  was  treated  by  allowing  all  par¬ 
ticles  to  simultaneously  expand  until  an  adequate  degree  of  overlap 
was  achieved.  Sunde  treated  the  structure  as  a  random  resistor 
network  and  use  Kirchhoff  s  current  laws  to  calculate  the  effec¬ 
tive  properties  of  the  microstructures.  For  validation,  results  were 
compared  to  independent  conductivity  measurements. 

Following  Sunde’s  efforts,  Martin  and  coworkers  tried  to 
improve  upon  the  accuracy  and  physical  relevance  of  these  types 
of  methods  through  the  use  of  a  discrete  element  method  (DEM), 
which  uses  Monte  Carlo  methods  to  generate  an  initial  structure 
followed  by  a  simulated  relaxation  process  to  represent  the  actual 
processing  steps  that  the  heterogeneous  SOFC  electrode  architec¬ 
ture  undergoes  [59,111-113].  This  approach,  first  used  to  describe 
particle  interactions  during  compaction  [59],  predicts  the  evolu¬ 
tion  of  particle  trajectories  according  to  Newton’s  laws  of  motion. 
Inter-particle  collisions  are  treated  during  relaxation.  Particles  are 
not  allowed  to  overlap  during  relaxation;  however,  they  are  sub¬ 
jected  to  a  volumetric  expansion  to  obtain  the  desired  contact  angle 
that  is  deemed  representative  of  the  sinter  process  in  a  manner  to 
that  of  Sunde.  A  recent  study  compared  DEM,  standard  Monte  Carlo, 
and  percolation  theory  predictions,  found  comparable  predictions 
of  coordination  number,  percolation  probability,  TPB  lengths,  and 
effective  conductivities  over  a  range  of  volumetric  compositions 
[114]. 

A  unique  stochastic  approach,  designed  for  the  reconstruction  of 
random  heterogeneous  and  multifunctional  media,  has  been  devel¬ 
oped  by  Torquato  and  coworkers  [115-119].  Torquato’s  approach 
assumes  that  real  structures  can  be  used  to  ascertain  correlation 
functions  that  represent  the  dispersion  and  interactions  of  the 
unique  phases  within  a  multifunctional  heterogeneous  material. 
If  a  correlation  function  adequately  describes  the  heterogeneous 
structure,  it  may  then  be  used  to  produce  an  infinite  number  of 
artificial  structures  that  are  statistically  representative  of  the  origi¬ 
nal.  This  approach  could  allow  an  infinite  number  of  permutations 
of  a  valid  microstructure  to  be  numerically  produced,  which  is  of 
great  use  when  trying  to  understanding  the  nature  and  role  of  the 
structure. 

As  a  part  these  efforts,  Torquato  has  shown  procedures  to  recon¬ 
struct  random,  multi-phase,  anisotropic  heterogeneous  media  from 
morphological  descriptions  of  a  base-line  system  [115].  This  is 
accomplished  using  an  n- point  correlation  function  to  provide  the 
degree  of  morphological  information  as  needed.  Reconstruction 
of  an  isotropic  ceramic-metallic  composite  (cermet)  was  demon¬ 
strated  by  using  orthogonal  sampling  with  a  two-point  correlation 
function  that  did  not  exhibit  short  range  order  effects.  However,  in 
addressing  uniqueness,  variability  between  reconstructions  with 
different  initial  guesses  is  noted  as  a  direct  result  of  finite  size 
effects  [117].  This  was  later  identified  as  a  more  significant  chal¬ 
lenge,  when  a  2-point  correlation  function  was  found  to  work 
well  for  reconstructing  heterogeneous  materials  with  a  single¬ 
scale  structure,  but  was  not  as  well  suited  to  a  multi-scale  random 


media  [118,119].  This  suggests  that  challenges  remain  for  using  this 
approach  with  broader  multifunctional  heterogeneous  materials, 
which  can  contain  hierarchical  structuring  with  the  length  scale  of 
features  spanning  several  orders  of  magnitude.  Despite  these  chal¬ 
lenges,  the  concept  of  statistical  correlation  to  actual  structures  is 
quite  interesting  and  unique. 

2.3.  Micro-  to  nano-scale  approaches 

The  consideration  of  discrete  heterogeneous  nano-  and  micro¬ 
structures  within  the  context  of  modeling  and  simulation  can 
provide  unique  details.  In  many  instances,  details  resolved  at 
nanometer  length  scales  (e.g.,  grains,  grain  boundaries,  inclusions, 
phase-interfaces,  contact/dihedral  angles,  and  pore  cross-sections) 
can  be  considered  while  also  maintaining  long-range  order  (i.e., 
properties  become  volume-independent)  that  is  often  found  in  the 
vicinity  of  several  to  tens  of  micrometers.  While  this  requires  the 
use  of  representative  volumes,  it  permits  the  unique  functional  pro¬ 
cesses,  and  their  coupling  to  the  discrete  heterogeneous  structure, 
to  be  captured.  The  physics  and  chemistry  must  be  understood  and 
defined  prior  to  simulation;  however,  the  influence  of  the  hetero¬ 
geneous  structure  on  these  processes  can  be  accounted  for  on  a 
non-empirical  basis  with  confidence  in  these  descriptions. 

There  have  been  numerous  groups  working  at  the  levels  of  the 
heterogeneous  structures  at  the  nanometer  to  micrometer  length 
scales.  Reifsnider  and  coworkers  examined  the  relative  contribu¬ 
tions  of  grain  boundary  diffusion  to  the  ionic  conductivity  of  ScSZ 
[17]  using  input  from  a  scanning  electron  microscopy  micrograph, 
where  elemental  phases  and  individual  grains  and  grain  bound¬ 
aries  could  clearly  be  distinguished.  This  study  was  performed  to 
highlight  the  importance  of  the  individual  grains  and  structure  at 
nanometer  length  scales  for  this  functional  heterogeneous  struc¬ 
ture.  Approximately  12%  difference  in  effective  conductivity  was 
identified  due  to  the  consideration  of  grain  boundary  diffusion 
mechanisms. 

The  study  by  Reifsnider  and  coworkers  highlighted  the  impor¬ 
tance  of  the  grain  structure  at  these  levels.  Efforts  by  Zhao  and 
Virkar  took  a  related  approach,  but  with  unique  considerations 
[37].  These  considerations  included  accounting  for  the  morphol¬ 
ogy  of  the  inter-particle  necks  and  the  contributions  of  discrete 
grain  boundaries  to  the  space  charge  regions  that  promote  defect 
stacking  and  can  assist  and/or  hinder  the  bulk  conductivity.  Zhao 
and  Virkar  considered  the  conductivity  of  porous  SDC  that  was 
made  via  standard  powder  compaction  method  and  using  NiO-SDC 
sintering.  The  NiO  was  subsequently  reduced  to  Ni  and  leached 
with  acid.  Porosity  of  the  samples  was  held  consistent;  however, 
the  acid-leached  SDC  sample  had  an  experimentally  measured  DC 
conductivity  that  was  up  two  orders  of  magnitude  larger  than 
the  compacted  sample.  Using  an  analytic  development,  these  dis¬ 
crepancies  were  attributed  to  the  inter-particle  neck  size,  which 
were  clearly  different  for  the  unique  methods.  In  addition  to  the 
inter-particle  neck  sizes,  space  charge  regions  contributions  were 
considered  for  both  dense  and  porous  samples.  Space  charge  contri¬ 
butions  were  identified  as  significant  in  porous  samples  with  small 
inter-particle  necks. 

The  contribution  of  space  charge  regions,  or  the  depletion  and/or 
stacking  of  electronic  and  ionic  defects  near  interfaces  in  solid 
state  devices,  is  a  topic  that  has  been  well  studied  due  its  impor¬ 
tance  to  nanostructured  systems  and  its  corresponding  interfaces 
[38,1 20-125].  Maier  et  al.  have  established  an  expertise  in  this  area, 
forming  both  analytic  and  numeric  approaches  to  describing  the 
effects  of  space  charge.  This  includes  aspects  of  the  all-important 
defect  chemistry  in  the  interfacial  region  [121-123]  and  the  associ¬ 
ated  thermodynamics  [  1 24].  These  approaches  have  been  extended 
into  chemical  transport  [38,121]  and  equivalent  circuit  models 
[125]  capable  of  describing  discrete  structures,  to  understand  the 
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geometric  morphological  affects  [38,120,121,123,124].  The  impor¬ 
tance  of  these  contributions,  as  recognized  by  Zhao  and  Virkar 
[37],  is  that  they  can  have  a  very  real  contribution  to  the  overall 
behavior  of  the  system  and  are  a  key  aspect  into  understanding 
this  behavior.  Further,  a  relatively  broad  class  of  test  materials 
coined  heterostructures  has  shown  anomalous  ionic  conductivi¬ 
ties  as  the  space  charge  regions  (i.e.,  Debye  length)  are  confined 
in  controlled  lattices,  which  have  been  validated  using  a  wealth 
of  Arrhenius  conductance  measurements  with  unique  structures 
[  1 20].  These  effects  become  more  pronounced  as  these  conductions 
zones  are  further  confined.  The  role  of  lattice  strain  in  the  vicin¬ 
ity  of  the  interface  may  also  be  questioned,  among  other  things. 
However,  the  measurement  of  drastic  increase  in  ionic  conduc¬ 
tivity  via  exploitation  of  these  principles  point  to  their  potential 
impact  [126].  Certainly,  these  principles  can  be  implemented  in 
stable,  functional  heterogeneous  materials  where  they  can  enhance 
conductive  transport,  including  at  lower  temperatures  which  offer 
a  number  of  cursory  benefits.  Such  an  undertaking  requires  both 
advances  in  scientific  understanding  of  the  physical  processes 
occurring  in  these  structures  to  relate  these  principles  into  rational 
model-based  design. 

LBM,  a  scheme  using  the  discretized  form  of  the  Boltzmann 
equation,  represents  a  simplified  kinetics  approach  to  describ¬ 
ing  transport  phenomena  in  a  discrete  heterogeneous  structure 
[62].  LBM  uses  a  velocity-space  representation  of  the  fluid,  using  a 
probability  distribution  function  (PDF)  to  provide  a  statistical  rep¬ 
resentation  of  the  fluid  and  its  interactions.  PDF  moment  closures 
are  used  to  recover  macroscopic  fluid  properties  including  con¬ 
centration  (or  density)  and  fluxes  (momentum,  energy).  Molecular 
interactions,  which  dictate  the  fluid’s  viscous  and/or  diffusive  pro¬ 
cesses,  are  treated  with  a  Bhatnagar-Gross-Krook  (BGK)  collision 
operator(s)  to  relax  the  system  versus  a  Maxwellian  distribu¬ 
tion.  Because  it  is  rooted  in  Boltzmann  statistics,  LBM  provides 
a  more  fundamental  description  of  molecular  interactions  than 
the  Navier-Stokes  equations,  which  can  be  exactly  recovered  [62]. 
Efforts  to  extend  the  LBM  formulation  through  the  transition  trans¬ 
port  regime  to  finite  Knudsen  numbers  have  also  been  undertaken 
[63-65].  Regular  meshes  at  rather  low  mesh  density  can  be  used 
with  considerable  fidelity,  which  is  a  substantial  benefit  for  han¬ 
dling  the  complex  features  observed  in  actual  structures  that  can 
make  traditional  meshing  schemes  nearly  intractable.  Transient 
relaxation  toward  equilibrium  occurs  as  an  explicit  time  march¬ 
ing  process.  The  PDF  at  each  lattice  point  is  only  dependent  upon  its 
nearest  neighbor  lattice  points,  which  are  updated  during  a  stream¬ 
ing  process  occurring  at  every  time  step.  Consequently,  the  explicit 
method  provides  nearly  linear  scalability  with  parallelization. 

LBM  has  been  used  to  perform  a  variety  of  numerical  simula¬ 
tions  within  the  SOFC  community.  Chiu  and  coworkers  have  used 
LBM  to  study  mass  and  charge  transport  in  the  heterogeneous 
structure  of  the  porous  Ni-YSZ  cermet  SOFC  anode  [64,127-135]. 
These  efforts  have  included  the  development  of  a  N-component 
multi-component  LBM  approach,  capable  of  simulating  contin¬ 
uum  [127]  and  non-continuum  [64]  diffusion  processes;  which 
were  subsequently  applied  to  2D  conceptual  electrode  structures  to 
ascertain  the  links  between  porosity,  tortuosity,  and  concentration 
polarizations  [128].  Initial  efforts  in  these  studies  were  validated 
using  independent  analytic  methods  and  experimentally  measured 
concentration  polarizations.  Asinari  et  al.,  used  LBM  with  granu¬ 
lometry  theory  based  scanning  electron  microscopy  micrographs 
of  actual  samples  to  simulate  transport  phenomena  in  a  recon¬ 
structed  3D  porous  Ni-YSZ  cermet  anode  to  calculate  localized 
fluxes  are  shown  and  macroscopic  tortuosity,  as  a  function  of  posi¬ 
tion  in  the  anode  [136].  The  efforts  of  Chiu  et  al.  and  Asinari  et  al. 
represent  some  of  the  first  efforts  to  utilize  reconstructions  of  actual 
heterogeneous  electrode  structures  to  perform  discrete  transport 
simulations. 


Mass  and  charge  transfer  were  also  explored  by  Chiu  and 
coworkers  in  actual  SOFC  anode  microstructures,  obtained  using 
tomographic  reconstruction  of  advanced  transmission  X-ray 
microscopy  data  and  supplemented  by  sphere-packing  gener¬ 
ated  structures  [130,131,133,134].  Nonlinearities  and  variations 
in  pore-scale  concentrations  and  potentials  were  observed  in 
these  studies,  resulting  from  the  complexity  of  the  actual  struc¬ 
ture.  The  power  of  these  methods  lie  in  their  capability  to  link 
morphology  to  their  contributions  to  transport  and  interfacial 
phenomena,  which  were  captured  when  an  explicit  quantita¬ 
tive  coupling  was  made  between  the  respective  grain  “particle 
sizes”  and  their  contributions  to  transport  losses  [133,134].  These 
details  were  captured  in  conjunction  with  quantitative  descrip¬ 
tions  of  the  size  distribution,  interfaces,  connectivity,  and  tortuous 
nature  of  the  morphology  to  demonstrate  the  connection  between 
these  morphological  descriptions  and  the  relative  performance 
[130,131,133,134,137].  Shikazono  et  al.  demonstrated  the  use  of 
FIB-SEM  reconstructed  RVEs  of  a  thin  porous  Ni-YSZ  cermet  anode 
to  incorporate  full  details  of  polarization,  including  electronic  and 
ionic  charge  transfer  processes  [138].  The  Butler-Volmer  equation 
was  used  to  represent  electrochemical  oxidation  and  the  electronic 
and  ionic  charge  transfer  processes  were  on  the  basis  of  gradients 
in  their  chemical  potential.  Following  validations  using  experimen¬ 
tal  measurements  of  button  cell  anode  polarizations  as  a  function 
of  current  density  with  different  inlet  H2  concentrations,  localized 
concentrations/potentials,  fluxes,  and  current  densities  were  iden¬ 
tified.  These  studies  demonstrated  the  considerable  nonlinearity 
of  the  overpotential— as  is  primarily  dictated  by  the  resistive  elec¬ 
trolytic  phase. 

As  demonstrated  by  the  Shikazono  et  al.,  one  of  the  many 
utilities  of  LBM  is  its  ability  to  capture  reactive  transport.  These 
processes  are  especially  important  in  the  SOFC,  because  they  couple 
the  discrete  functional  processes  through  the  heterogeneous  struc¬ 
tural  interfaces.  Additional  studies  have  been  performed  by  Chiu 
et  al.  to  explicitly  link  elementary  kinetics  mechanisms  describing 
detailed  electrochemical  oxidation  kinetics  [129,132]  and  detailed 
heterogeneous  internal  methane  reformation  kinetics  [129,135]  to 
the  transport  processes  to  better  understand  the  effect  of  this  cou¬ 
pling.  These  approaches  showed  the  sensitivity  of  the  coupling,  role 
and  importance  of  the  discrete  microstructure,  and  unique  kinetic 
regimes  of  reactive  transport  in  these  systems.  Such  trends  are 
characteristic  of  a  multi-scale  coupled  system. 

Phase  field  modeling  is  also  an  important  tool  that  has  shown 
considerable  insights  into  the  aging  and  the  evolution  of  the  SOFC 
microstructures,  which  occur  over  broad  time  scales.  These  mod¬ 
els  can  describe  microstructural  evolution  due  to  processes  like 
Ostwald  ripening  [139],  precipitation  reactions  [140],  and  heat 
treatment  and  thermally  induced  coarsening  [141 ,142],  among  oth¬ 
ers.  Empirical  inputs  and/or  calibration  are  typically  required  as 
not  all  of  the  degradation  mechanisms,  rates  and  processes  are 
quantitatively  known  due  to  the  time  scales  involved. 

This  is  consistent  with  a  work  by  Kim  et  al.,  who  present  a  phase- 
field  modeling  approach  for  the  SOFC  anode  [143].  Aging  of  the 
SOFC  anode  microstructure  was  described  using  the  Cahn-Hilliard 
equation,  which  is  calibrated  with  experimental  input  to  kinetic 
parameters.  The  Cahn-Hilliard  equation  in  just  one  approach,  as 
analytic,  boundary  integral  and  general  phase  field  approaches  are 
also  available.  The  Cahn-Hilliard  equation  is  used  to  describe  con¬ 
tinuity  of  material  flux  due  to  differences  in  free  energy  across  an 
interface.  Sinusoidal  representations  of  phase  interfaces  in  equilib¬ 
rium  were  used  and  curvature  corresponding  to  particle  sizes  pro¬ 
vides  additional  free  energy  contributions.  Cell  performance  was 
described  with  homogenized  empirical  parameters  in  representa¬ 
tive  volume  elements.  A  series  of  representative  volumes  were  used 
to  discretize  the  anode  thickness  and  the  associated  mass  transfer, 
charge  transfer,  and  phase  field  representations  were  solved  using 
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a  finite  difference  method.  Phase-field  variables  were  assigned  to 
each  RVE  and  changes  in  particle  diameter  and  TPB  length  were 
studied  over  1 00,000  h  of  operation.  Changes  in  Ni  particle  diame¬ 
ter  matched  experimentally  measured  values  obtained  over  3000  h 
of  cell  operation.  The  experimental  measurements  were  taken  from 
the  literature,  which  consisted  of  characterization  made  using  line 
intercept  methods  on  micrographs  of  sectioned  samples. 

In  a  work  demonstrating  the  capabilities  enabled  by  having 
physical  representation  for  the  3D  microstructure,  Chen  et  al. 
used  FIB-SEM  reconstructions  of  the  Ni-YSZ  cermet  SOFC  anode 
as  an  input  to  a  phase  field  model  [144].  For  this  effort,  the 
group  studied  the  coarsening  of  Ni  in  this  system  with  aging 
and  more  specifically  the  implications  of  assumptions  regarding 
the  mobility  of  the  secondary  phases  and  correlate  phase  field 
results  to  physical  measurements  of  the  microstructure,  such 
as  changes  in  electrochemically  active  TPB  length,  inter-phase 
contact  area,  and  tortuosity.  While  this  particular  study  did  not 
perform  detailed  validations,  these  are  properties  that  have  been 
measured  in  the  SOFC  by  others  using  FIB-SEM  and  X-ray  based 
tomographic  methods  [131,133,134,137,145-148].  Interestingly, 
Chen  et  al.  [144]  note  the  importance  of  the  contact  angles  on 
the  aging  behavior,  with  a  reduction  in  the  Ni-YSZ  contact  angle 
decreasing  reconfiguration  within  the  microstructure.  This  is  not 
unexpected,  as  the  contact  angle  dictates  the  force  of  adhesion 
between  the  phases  and  the  energy  that  is  ultimately  required  for 
reconfiguration;  however,  the  ability  to  correlate  and  map  these 
effects  between  the  theoretical  model  and  physical  observations 
in  the  actual  system  promise  to  greatly  improve  the  fidelity  and 
capability  of  these  phase  field  methods  and  provide  access  to 
suggest  materials  and  design  changes. 

Besides  aging  and  microstructural  evolution,  there  is  merit  to 
considering  the  mechanical  aspects  of  the  grain-level  structure. 
There  are  some  considerable  challenges  to  such  an  effort,  such 
as  applying  boundary  conditions,  accounting  for  residual  stresses, 
and  resolving  an  accurate  and  tractable  thermal  model  if  aspects 
such  as  thermal  expansion  want  to  be  considered.  With  a  few 
exceptions,  these  types  of  challenges  have  significantly  limited  the 
work  done  in  this  area.  Micromechanical  models  describing  frac¬ 
ture  in  heterogeneous  and  grain-level  systems  are  available  and 
typically  limited  to  cohesive  fracture  in  brittle  materials  [9,149]. 
These  approaches  again  use  RVEs,  where  actual  grain  morpholo¬ 
gies  are  required.  Within  the  RVE,  finite  element  descriptions  of  an 
anisotropic  elastic  media  are  used  with  cohesive  interface  elements 
along  the  crack  front.  As  noted  in  these  studies,  the  difficulty  lies  in 
validating  the  models  and  approaches  as  they  inherently  occur  at 
relatively  random  time  and  locations  in  the  actual  system.  Aspects, 
such  as  crack  patterns  and  velocity  histories  can  provide  insight,  but 
the  challenge  remains  formidable.  Still,  the  benefit  of  such  models 
is  being  able  to  predict  constitutive  material/structural  properties 
both  with  and  without  defects. 

Possibly  a  more  tractable  example  and  demonstration  of 
changes  in  the  mechanical  behavior  of  the  heterogeneous  electrode 
structures  has  been  demonstrated  by  Liu  et  al.,  where  the  effect  of 
phosphorous  poisoning  on  the  structural  integrity  of  the  cell  was 
explored  [54].  Finite  elements  are  used  to  discretize  the  linear  elas¬ 
tic  and  constitutive  equations  within  actual  2D  SEM  micrographs 
of  sectioned  porous  Ni-YSZ  cermet  anodes  that  have  been  exposed 
to  different  levels  of  phosphorus.  For  these  efforts,  electron  disper¬ 
sive  spectroscopy  (EDS)  was  used  to  identify  Ni  and  YSZ,  as  well  as 
the  conversion  of  Ni  to  NixPy  phases.  Using  the  finite  element  anal¬ 
ysis  on  RVEs  of  the  structure,  von  Misses  stresses  in  the  Ni  and  YSZ 
phases  of  the  structure  were  examined,  along  with  their  respective 
evolution.  These  stresses  were  the  result  of  volumetric  changes  of 
Ni  following  phosphorous  poisoning.  Due  to  the  strong  interactions 
of  Ni  and  P,  contamination  followed  a  wave-front  type  propagation 
with  exposure  and  a  unique  stress  evolution  with  time. 


2.4.  Molecular  simulations 

Molecular  dynamics  (MD)  can  provide  a  great  deal  of  informa¬ 
tion  about  the  rates  of  processes  and  the  relaxation  of  materials. 
From  a  theoretical  standpoint  MD  is  rather  straight-forward;  New¬ 
ton’s  laws  of  motion  are  applied  to  particles  in  the  system,  which 
correspond  to  atoms  and  molecules.  Inter-particle  interactions  can 
take  the  form  of  electrostatic,  Van  der  Wall’s,  and  Lennard-Jones 
potentials,  as  well  as  the  various  modes  of  atomic  and  molec¬ 
ular  motion.  Many  implementations  use  force  fields  to  describe 
inter-atomic  and  molecular  potentials  that  are  either  long-range 
or  involve  more  complex  processes  that  are  dependent  upon  the 
atom/molecule’s  electronic  structure.  In  terms  of  using  MD  meth¬ 
ods  to  study  functional  heterogeneous  structures,  MD  methods  are 
quite  suited  to  finite  periodic  regions  descriptive  of  an  interface  or 
several  unit  cells  of  structure  at  this  time  due  to  their  computational 
costs. 

Among  the  more  common  approaches  used  within  the  SOFC 
community  is  a  kinetic  Monte  Carlo  (KMC)  method.  Prinz’s  group 
has  made  use  of  KMC  for  studying  the  mechanisms  and  proper¬ 
ties  of  electrolyte  materials  [150,151].  These  KMC  studies  require 
appropriate  inter-atomic  and  inter-molecular  descriptions,  as  dis¬ 
cussed  in  the  following  section.  In  these  studies,  KMC  has  been 
used  to  develop  an  electrochemical  impedance  type  simula¬ 
tion  of  the  material’s  response  and  specifically  details  in  regard 
to  the  oxygen  diffusion  at  the  space  charge  double  layer  at 
the  electrode-electrolyte  interfaces  and  the  role  of  the  elec¬ 
trolyte  thickness  [150].  This  enabled  experimental  validation  to 
impedance  spectra.  Information  about  the  attributed  energy  barri¬ 
ers  for  electrolyte  materials  were  required  from  density  functional 
theory  simulations.  Using  alternating  potentials  to  represent  the 
impedance  response,  low  frequency  oxide  migration  was  found  to 
be  fast  in  the  electrolyte  and  that  there  was  an  accumulation  of 
oxide  at  the  electrode-electrolyte  interfaces  forming  a  double  layer 
and  space  charge  region.  Thicker  electrolytes  were  found  to  require 
longer  times  to  obtain  the  same  double-layer  characteristics.  Only 
local  ionic  movement  was  recognized  at  higher  frequencies.  Double 
layer  capacitances  extracted  from  the  simulations  were  found  to  be 
considerably  larger  than  geometric  capacitance  contributions. 

In  a  subsequent  study  using  simulation  methods  to  fur¬ 
ther  understand  impedance  spectra,  Prinz  and  coworkers  used 
KMC  to  model  vacancy  migration  in  the  solid  electrolyte  to 
show  that  the  fluctuation  dissipation  theorem  holds  across  all 
frequencies  [151].  This  study  examined  the  effects  of  dopant 
concentrations  in  the  oxide  electrolyte.  It  was  found  that  at  infi¬ 
nite  frequency,  impedance  decreases  with  increased  doping.  And 
at  zero-frequency,  dopant-vacancy  interactions  are  dominant. 
Approaches  to  extend  these  efforts  for  the  consideration  of  gen¬ 
eralized  surfaces  and  grain  boundary  contributions  are  discussed. 

Goddard  III  and  coworkers  have  performed  MD  based  studies 
of  oxygen  ion  transport  in  YSZ  [152].  Reactive  force  fields  have 
been  developed  and  optimized,  which  are  used  to  study  oxygen 
ion  diffusion  in  YSZ  as  a  function  of  temperature.  The  results  were 
validated  with  experimentally  measured  Arrhenius  data,  with  high 
fidelity.  This  group  notes  that  they  are  trying  to  develop  an  accurate 
description  of  oxide  transport  in  YSZ  so  that  force  fields  for  various 
catalysts  can  be  used  to  enable  simulation  of  critical  processes  of 
these  materials  and  interfaces  and  subsequent  optimization  at  the 
interface-level. 

Moving  to  the  examination  of  the  oxygen  reduction  reaction 
(ORR)  in  the  SOFC  cathode,  Dunlap  and  corkers  use  KMC  to  study 
the  mechanisms  and  rates  [153-155].  First,  KMC  was  used  to  study 
the  LSM-YSZ  interface  and  its  affect  on  the  ORR  [153].  A  sensitiv¬ 
ity  analysis  using  the  ionic  current  density  was  used  for  this  effort, 
and  was  found  to  be  most  sensitive  to  variations  of  the  O2-  incor¬ 
poration  rate  into  the  YSZ.  This  is  a  result  of  a  large  oxygen  vacancy 
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diffusion  barrier.  Oxygen  segregation-induced  concentration  gra¬ 
dients  at  the  double  layer  were  recognized  as  having  key  intrinsic 
properties  that  affect  performance.  In  a  subsequent  study,  Dunlap 
and  coworkers  also  examine  both  the  ORR  with  a  Pt-YSZ  cathode 
and  the  hydrogen  oxidation  reaction  (HOR)  of  the  Ni-YSZ  anode 
[154].  This  study  identified  oxygen  incorporation  into  the  YSZ  elec¬ 
trolyte  as  rate  limiting.  Steps  also  of  importance  but  will  a  smaller 
impact  on  the  overall  rate  were  identified  as  being  the  oxide  ion 
transport  in  YSZ  and  oxygen  spillover  onto  Ni  for  the  HOR. 

An  interesting  possibility  for  molecular  level  simulations  that 
has  not  been  extensively  used  in  the  context  of  heterogeneous  SOFC 
is  the  use  of  such  methods  for  examining  the  mechanical  aspects 
of  the  constituent  materials  and  interfaces.  Sato  et  al.  reported  the 
use  of  MD,  along  with  experimental  efforts,  to  study  the  effects  of 
compressive  and  tensile  strain  on  the  conductivity  of  Y203  doped 
Ce02,  where  different  amounts  of  dopant  were  considered  [156]. 
These  studies  identified  peak  conductivity  at  a  dopant  concentra¬ 
tion  of  approximately  20%  Y203,  which  was  in  agreement  for  both 
simulation  and  experimental  efforts.  Elastic  moduli  were  also  in 
agreement  and  scaled  inversely  with  conductivity. 

Crack  formation  and  propagation  dynamics  are  topics  that  have 
been  extensively  studied  at  the  atomic  and  molecular  levels;  how¬ 
ever,  do  not  appear  to  have  ventured  very  far  into  the  SOFC 
community  [7,8].  The  complex  interactions  of  time/length  scales 
and  materials/structures  involved  presents  a  computational  chal¬ 
lenge  for  addressing  these  issues  in  the  SOFC.  As  they  become  more 
tractable,  such  studies  can  provide  details  regarding  the  dynam¬ 
ics,  pathways,  and  rates  of  fractures  in  complex  systems  and  have 
specifically  been  applied  to  both  brittle  and  ductile  materials  [8]. 
These  approaches  are  of  particular  use  for  examination  of  highly 
strained  and  broken  bonds,  where  quantum  approaches  are  neces¬ 
sary  for  higher  order  accuracy  of  these  processes  [7]. 

2.5.  Electronic  structure  and  ab  initio  methods 

The  electronic  structure  of  the  atoms  and  molecules  of  the  par¬ 
ticipating  materials  and  interfaces  can  provide  insight  into  the 
properties  and  nature  of  the  system.  Many  of  the  configurations, 
properties,  and  mechanistic  pathways  may  be  addressed  by  cal¬ 
culation  of  the  electronic  structure  of  ensembles  of  atoms  and 
molecules  and  the  minimization  of  the  energy  associated  with 
the  structure/ensemble.  As  with  molecular  simulations,  electronic 
structure  calculations  are  typically  too  expensive  to  consider  a 
large  volume.  They  are  better  suited  to  extracting  properties  and 
details  regarding  the  constituent  materials,  interfaces,  configura¬ 
tions,  energetic,  and  processes/pathways,  which  can  be  described 
in  a  periodic  system.  These  facets  are  all  important  to  interfacial 
phenomena  like  transport  mechanisms/rates,  reaction  pathways, 
preferential  adsorption  sites/configuration,  activation  energies, 
mechanical  properties,  chemical  stability,  etc.,  which  can  calculated 
using  minimum  energy  pathways  [157]. 

Traditional  DFT  methods  use  a  ground-state  energy  density  that 
satisfies  Schrodinger’s  equation  to  predict  properties  via  treatment 
of  the  electronic  density  and  ground  state  energy  [55,157].  This 
may  be  achieved  using  a  Kohn-Sham  orbital,  which  express  the 
ground-state  energy  as  a  function  of  the  kinetic  energy,  exchange 
correlation  energy,  and  interactions  between  the  nuclei,  electrons, 
and  themselves.  Tight  binding  methods  may  be  used  to  separate 
contributions  between  core  and  valence  electrons,  the  latter  of 
which  dominates  covalently  bonded  systems  which  traditional  DFT 
methods  excel  with  [34].  The  Kohn-Sham  orbital  requires  a  suitable 
exchange  correlation  energy  term.  Local  density  (LDA)  and  gen¬ 
eralized  gradient  approximations  (GGA)  are  often  used.  The  GGA 
is  noted  for  its  balance  of  cost  and  providing  qualitatively  correct 
chemical  properties,  making  it  a  method  of  choice  for  studying 
topics  like  heterogeneous  catalysis  [1].  Suitable,  functional  forms 


are  required  for  the  generalized  gradient  approximation,  such  as 
Perdew-Burke-Ernzerhof  (PBE)  and  Perdew-Wang  91  (PW91). 

In  the  SOFC  community,  a  number  of  research  groups  have 
been  using  DFT-based  methods  to  understand  material,  interface, 
and  nanostructure  properties  as  well  as  physical/chemical  path¬ 
ways.  In  terms  of  the  conductivity  of  the  common  YSZ  electrolyte 
material,  Pornprasertsuk  et  al.  used  a  DFT  method  using  LDA  and 
gradient  correction  with  kinetic  Monte  Carlo  to  study  the  oxygen 
vacancy  diffusion  mechanism  [158].  The  effect  of  dopant  concen¬ 
tration  on  conductivity  was  examined,  where  DFT  provided  the 
energy  barriers  along  the  minimum  energy  path  for  a  vacancy  to 
migrate/diffuse.  Traditional  material  defect  theory  would  suggest 
that  conductivity  should  continue  to  increase  with  dopant  concen¬ 
tration;  however,  it  has  long  been  recognized  that  the  conductivity 
no  longer  increases  once  approximately  8mol%  Y203  is  reached. 
Pornprasertsuk  et  al.  identified  the  migration  energy  needed  to 
traverse  two  adjacent  tetrahedral  containing  the  Zr-Y  and  Y-Y 
common  edge  as  the  root  of  this  limitation  [158].  Experimental  and 
simulated  conductivity  versus  dopant  concentration  trends  were  in 
strong  accord. 

Adsorption,  dissociation,  and  pathways  associated  with  cat¬ 
alytic  and  electrocatalytic  processes  supported  by  the  heteroge¬ 
neous  electrode  structures  were  explored  by  Rosmeisl  and  Bessler, 
whom  used  plane-wave  DFT  with  GGA  to  calculate  the  stability  of 
H,  0,  and  OH  radicals  on  Mn,  Fe,  Co,  Ni,  Cu,  Ru,  Rh,  Pd,  Ag,  Pt,  and  Au 
[  1 59].  These  materials  were  treated  as  candidate  anode  electrocata¬ 
lyst  materials  and  were  interfaced  with  YSZ  so  that  electrochemical 
oxidation  pathways  could  be  examined  via  energy  minimization. 
Oxygen  spillover  was  identified  as  a  dominant  pathway,  with  the 
electrocatalytic  activity  linked  to  the  binding  affinities  of  the  reac¬ 
tive  intermediates.  The  electrocatalytic  activity,  defined  as  the 
change  in  free  energy  for  the  rate  determining  step,  was  found  to 
scale  with  the  experimentally  determined  binding  energy  for  the 
various  catalysts  with  Ni  sitting  atop  the  volcano  plot  [159]. 

Like  Rosmeisl  and  Bessler,  Shishkin  and  Ziegler  [160]  and  Ingram 
and  Linic  [161]  used  similar  approaches  to  study  the  activity  of 
various  materials  for  direct  hydrocarbon  oxidation.  Shishkin  and 
Ziegler  used  a  DFT  algorithm  employing  a  PBE  based  exchange- 
correlation  functional  using  spin-polarized  calculations  and  a 
projected  augmented  wave  method  to  study  the  oxidation  of  H2, 
CH4,  and  CO  at  the  Ni/YSZ  interface  (i.e.,  TPB)  with  a  focus  on  the 
charge  transfer  processes.  In  this  effort,  Shishkin  and  Ziegler  noted 
that  H2  and  CH4  have  relatively  large  adsorption  activation  bar¬ 
riers,  whereas  CO  is  favorable  but  endothermic.  Oxygen  spillover 
from  YSZ  is  suggested  as  being  dominant,  but  it  is  further  suggested 
that  some  H-spillover  from  Ni  to  YSZ  may  compliment  this  process. 
Further,  it  is  observed  that  H  can  oxidize  an  O-enriched  surface, 
but  that  O-enriched  YSZ  is  less  active  near  Ni  than  an  infinite  YSZ 
surface  [160].  Ingram  and  Linic  used  a  plane  wave  basis  set  based 
DFT  with  electron  exchange  and  correlation  being  described  by  the 
PW91  function  with  a  GGA  to  study  various  metal  catalysts,  and 
how  they  behave  in  the  presence  of  hydrocarbons  [161  ].  They  note 
that  noble  metals  are  relatively  inert  to  coking  because  they  do  not 
activate  C-H  bonds  and  that  active  metals,  such  as  W,  Fe,  and  Mo 
are  poisoned  by  oxygen  and  carbon.  This  implies  that  metals  like 
Co,  Ru,  Ir,  Rh,  and  Ni  offer  the  optimal  performance  because  they 
can  activate  CH4  for  charge  transfer  while  not  necessarily  oxidiz¬ 
ing.  These  metals  still  suffer  from  C-poisoning,  depending  upon  the 
C/O  ratio  due  to  the  formation  of  C-0  rather  than  C-C  bonds  [161  ]. 

Besides  the  anodic  processes,  the  oxygen  reduction  reaction  in 
the  cathode  has  also  extensively  been  studied,  as  evident  by  the 
review  by  Choi  et  al.,  who  highlight  the  challenges  of  the  surface 
and  interface  models  for  the  junction  of  the  electrode  and  elec¬ 
trolytic  materials  [157].  In  this  review  one  of  the  topics  examined 
was  the  adsorption  sites  and  configurations  for  perovskite  cath¬ 
ode  materials,  as  well  as  the  oxygen  incorporation  step.  A  study 
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by  Mastrikov  et  al.,  whom  used  DFT  calculations  with  exchange 
correlations  using  a  GGA  with  a  PW91  function  and  plane  wave 
basis  set,  studied  the  mechanism  of  oxygen  incorporation  into 
(La!_xSrx)Mn03_5  [162].  As  a  part  of  this  effort,  Mn02  (001)  was 
identified  as  the  most  stable  surface  termination  for  fuel  cell  con¬ 
ditions.  Mechanistically,  the  rate  determining  step  was  identified 
as  the  surface  oxygen  vacancy  encountering  an  adsorbed  O-  via 
activated  surface  diffusion  of  the  vacancy.  This  points  toward  the 
need  for  a  higher  surface  vacancy  concentration  and  mobility.  As 
a  part  of  these  efforts,  oxygen  adsorbates  were  found  to  be  rather 
low,  which  was  explained  by  negative  adsorption  entropy  and  elec¬ 
trostatic  repulsion,  despite  the  exothermic  formation  of  02_,  022-, 
and  O-  on  Mn. 

Another  prominent  concern  with  the  materials  used  in  the 
SOFC’s  heterogeneous  materials  is  their  stability  in  the  presence 
of  poisons  and  contaminants.  In  an  attempt  to  address  the  stability 
and  degradation  mechanisms  of  Ni  in  the  presence  of  carbonaceous 
fuels,  Nikolla  et  al.  used  DFT  with  the  exchange  correlation  formu¬ 
lated  using  a  GGA  with  a  PW91  function  set  [163].  They  identified 
that  alloying  the  Ni  surface  can  promote  the  oxidation  of  C  atoms, 
rather  than  allowing  the  formation  of  strong  C-C  bonds,  which  have 
a  lower  driving  force  and  thus  permit  the  nucleation  of  C  atoms  on 
low-coordinated  Ni  sites.  In  conjunction  with  experimental  efforts, 
a  Sn/Ni  alloy  was  shown  to  improve  stability  to  carbon  and  adsorp¬ 
tion  relative  to  Ni  (1 1 1 ),  in  the  presence  of  methane,  propane,  and 
isooctane  with  a  steam  to  carbon  ratio  of  1 .5.  These  demonstrations 
experimentally  exhibited  hundreds  of  hours  of  stability.  Galea  et  al. 
performed  similar  simulations  on  Ni  and  Cu  (1 1 1)  and  (211),  as 
well  as  Cu-Ni  and  Cu-Co  alloys  in  the  presence  of  CFI4  [164].  They 
found  that  C  adsorbs  strongly  onto  the  Ni  (2 1 1 )  surface  and  grows 
graphitic  carbon  over  the  terrace;  however,  Cu  was  found  to  have 
a  very  high  thermodynamic  and  kinetic  barrier  to  CFI4  dissociation. 
This  dissociation  barrier  restricts  breakdown  and  provides  oppor¬ 
tunities  for  direct  electrochemical  oxidation.  The  Cu-Ni  and  Cu-Co 
alloys  showed  that  Ni  and  Co  has  little  effect  on  surface  Cu  and  that 
stability  is  limited  by  the  Cu  enrichment  of  the  alloy  surface. 

Besides  carbon,  sulfur  is  also  a  particularly  challenging  contami¬ 
nant  for  the  SOFC  anode.  Wang  and  Liu  used  DFT  with  the  exchange 
correlation  treated  using  a  GGA  with  the  PBE  function  to  study  the 
regeneration  of  a  Ni  catalyst  surface  with  02  and  H20  after  the  cat¬ 
alyst  was  exposed  to  S  [41  ].  The  study  found  that  both  02  and  H20 
were  practical  for  sulfur  removal;  however,  FI20  was  suggested  as 
the  better  choice  because  of  its  ability  to  remove  sulfur  over  a  wider 
range  of  pressures  without  oxidizing  the  Ni.  The  study  considered 
the  removal  and  adsorption  of  sulfur,  as  well  as  the  oxidation  of  the 
Ni  surface  to  form  NiO. 

Reduction  and  re-oxidation  (redox)  instability  in  the  SOFC  has 
been  widely  documented  as  leading  to  irreversible  degradation, 
typically  through  mechanical  changes  in  the  structure.  Coarsen¬ 
ing  and  irreversible  changes  have  been  suggested  as  the  culprit; 
however,  details  of  the  exact  mechanisms  remain  elusive.  In  a 
study  by  Jeangros  et  al.,  DFT  using  GGA  with  PBE  exchange  and 
correlation  has  been  used  to  confirm  in  situ  high  resolution  environ¬ 
mental  transmission  electron  microscope  (FIR-TEM)  reduction  and 
re-oxidation  (redox)  cycle  studies  of  the  Ni-YSZ  cermet  SOFC  anode 
[26].  This  investigation  linked  these  experiments  to  DFT  based  stud¬ 
ies  of  the  pathways  for  the  reduction  and  re-oxidation  of  the  Ni/NiO 
in  the  Ni-YSZ  cermet  anode.  For  this  effort,  vacancies  near  the 
Ni-YSZ  interface  were  introduced  with  the  nudged  elastic  band 
(NEB)  approach.  Experimentally,  reduction  showed  the  transfor¬ 
mation  of  NiO  grains  to  metallic  Ni,  leaving  a  nanoporous  structure. 
The  formation  of  additional  Ni  nanoparticles  was  also  observed. 
Jeangros  et  al.  suggested  that  the  reduction  process  initiated  with 
the  reduction  of  NiO  at  the  NiO/YSZ  interface  followed  by  that  of  the 
free  NiO  surfaces.  The  NEB  calculations  confirmed  this  behavior  as 
it  showed  that  the  relative  defect  energy  reduced  when  the  oxygen 


was  incorporated  into  the  YSZ,  which  was  related  to  relaxation  in 
strain  with  the  vacancy  in  the  NiO.  This  led  to  a  two  step  pathway 
being  proposed,  in  which  oxygen  ions  are  transferred  from  NiO  to 
YSZ,  which  subsequently  desorbs  as  water.  The  first  involved  the 
transfer  of  the  oxygen  to  the  YSZ  at  the  NiO/YSZ  interface,  followed 
by  the  reduction  of  the  oxygen  in  the  YSZ  with  hydrogen,  forming 
water  that  subsequently  desorbs.  The  authors  noted  that  a  comple¬ 
mentary  mechanism  of  hydrogen  adsorbing  onto  the  Ni  adjacent  to 
oxygen  vacancies  may  allow  water  to  desorb  directly  at  the  NiO/Ni 
interface  and  allow  NiO  reduction  to  continue  toward  the  center  of 
the  grain,  which  is  consistent  with  the  direct  reduction  of  the  free 
Ni/NiO  surface  [26]. 

While  GGA  based  DFT  approaches  have  been  mainstays,  DFT  +  U 
and  similar  hybrid  methods  are  becoming  more  prominent  because 
they  are  becoming  more  robust  and  permit  additional  relevant 
processes  and  materials  to  be  studied  [1,34].  Huang  and  Carter 
note  that  hybrid  schemes  like  DFT  +  U  provide  less  expensive  cal¬ 
culations  for  materials  involving  mixed  electron  distributions  and 
highly  correlated  materials  [1,34].  This  specifically  relates  to  struc¬ 
tures  and  materials  incorporating  transition  metal  oxides  and 
sulfides,  which  are  prevalent  in  the  SOFC  materials  and  processes. 
The  improvement  of  DFT  +  U  methods  results  from  its  exact  incor¬ 
poration  of  the  exchange  process. 

An  example  of  the  use  of  the  DFT  +  U  framework  is  a  work  by 
Chen  et  al.,  whom  used  the  method  to  study  the  sulfur  tolerance 
of  Ceria  (Ce02)  [165].  The  study  required  DFT  +  U  to  appropriately 
correct  for  the  on-site  coulomb  repulsion  of  the  Ce4f  state  local¬ 
izations.  The  general  calculations  used  GGA  based  DFT  with  PW91 
exchange-correlation  functional.  This  study  examined  the  interac¬ 
tions  of  hydrogen  sulfide  (H2S)  with  Ce02  (111),  with  H2S,  HS,  and 
S  intermediates,  which  were  produced  as  a  result  of  low  dehydro¬ 
genation  barriers.  Adsorption  sites  were  explored,  finding  a  variety 
of  bridging  sites.  Dehydrogenated  surface  adsorbates  were  further 
shown  to  react  with  available  oxygen  anions  under  operational  con¬ 
ditions,  with  an  S02  forming  pathway  being  the  most  likely  [165]. 
While  not  directly  compared,  surface  vibrational  spectra  were  cal¬ 
culated  as  a  part  of  this  effort  so  that  they  could  be  compared  to 
experimental  measurements  in  future  studies. 

Another  DFT  +  U  study  was  reported  by  Shishkin  and  Ziegler 
to  examine  the  chemical  properties  of  Ni/Ceria  (Ce02)  and 
Ni/Ce02/YSZ  SOFC  anodes  [166].  DFT  +  U  was  required  because  of 
the  coulomb  repulsion  of  the  Ce4f  electrons.  For  this  effort,  charge 
localization  resulting  from  vacancy  formation  was  examined,  as 
they  relate  to  the  adsorption  and  reaction  of  H2  and  CH4  fuels  on  the 
ceria  surface.  The  adsorption  of  H2  and  CH4  onto  the  Ni/Ce02  sur¬ 
face  approaches  those  of  pure  Ce02,  suggesting  similar  mechanisms 
and  rates.  Both  hydrogen  and  oxygen  spillover  pathways  were 
examined,  where  oxygen  spillover  onto  Ni  was  both  significant 
and  suggestive  of  the  higher  coking  tolerance.  During  oxidation, 
a  ceria  surface  vacancy  formation  resulted  in  the  transfer  of  a  small 
amount  of  charge  to  the  Ni.  This  is  different  than  Ni/YSZ  where  the 
valance  band  structure  dictates  that  it  is  energetically  favorable  to 
transfer  the  electronic  charge  to  the  Ni.  Still,  the  more  prominent 
pathway  was  the  oxidation  of  the  fuel  at  the  Ce02  surface  [166]. 

Recently,  Carter  et  al.  applied  DFT  +  U  methods  to  study  oxide 
ion  conductivity  in  perovskite-type  transition  metal  oxides 
LaM03  (M  =  Cr,  Mn,  Fe,  Co)  for  SOFC  cathode  materials  [168] 
and  Sr2FeMo06,  mixed  ionic  and  electronic  conductors  proposed 
as  electrodes  for  SOFC  applications  [169].  Bulk  oxygen  vacancy 
formation,  which  is  a  key  factor  for  the  oxide  ion  bulk  diffusion 
process,  was  carefully  examined  for  various  material  composi¬ 
tions.  Important  design  principles  were  derived  based  on  easily 
measurable  or  computable  properties. 

Several  additional  methods,  which  have  not  been  extensively 
used  in  the  SOFC  community,  include  quantum  Monte  Carlo  (QMC) 
and  AIMD  methods.  QMC  is  very  accurate  for  highly  correlated 
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systems,  but  remains  expensive  and  difficult  to  implement  [1,34]. 
Computational  expense  also  limits  AIMD  based  methods,  which 
link  molecular  mechanics  and  electronic  structure  calculations 
[1,34,55].  AIMD  may  be  considered  a  part  of  a  broader  class  of 
embedded  cluster  methods,  which  try  to  incorporate  longer  range 
effects  by  coupling  the  domain  with  more  general  theory  via  care¬ 
fully  selected  partitioning.  However,  this  becomes  difficult  as  size 
effects  become  important.  Certainly,  AIMD  has  had  limited  appli¬ 
cability  in  the  SOFC  community  to  date  because  of  the  time  and 
length  scales  required  to  perform  analysis  at  SOFC  conditions  with 
consistent  calculations  of  the  electronic  structure. 

3.  Conclusions  and  outlook 

Over  the  course  of  this  review,  a  variety  of  studies  that 
have  explored  unique  aspects  and  scales  of  the  SOFC  system, 
operation,  structure,  chemical  and  physical  processes,  and  degrada¬ 
tion/evolution  have  been  highlighted.  Many  of  these  studies,  which 
are  a  sampling  of  those  within  the  community,  have  taken  unique 
aspects  to  understand  processes  in  the  SOFC.  Individually,  each  of 
the  methods  presented  is  capable  of  providing  valuable  information 
regarding  the  properties,  performance,  structure,  and  interactions 
within  the  heterogeneous  structures  comprising  the  electrode 
structures  of  the  SOFC,  which  should  improve  with  advancements 
in  computational  facilities  and  theory/algorithms. 

These  individual  incremental  advances  hold  tremendous  value; 
however,  more  substantial  advances  may  come  from  the  cou¬ 
pling  of  these  approaches  and  scales  to  ascertain  more  of  the 
multi-scaled  nature  of  the  heterogeneous  SOFC  materials.  Incor¬ 
porating  the  strengths  of  these  individual  methods,  time  and 
length  scale  coupling,  along  with  those  of  the  unique  physical, 
chemical  and  mechanical  processes  that  constitutes  the  materi¬ 
als’  behavior  should  be  able  to  be  understood  at  unprecedented 
levels.  This  includes  predictive  capabilities  for  never-before  manu¬ 
factured  materials  [1  ].  Whether  automated  or  with  explicit  human 
interaction,  this  coupling  should  allow  inter-dependencies  of  the 
numerous  variables  involved  in  the  properties  and  function  of  the 
heterogeneous  material  system  to  be  systematically  understood.  A 
key  aspect  to  realizing  such  a  robust  level  of  modeling  and  simula¬ 
tion  is  the  complementary  development  of  improved  experimental 
capabilities;  specifically  those  that  can  provide  new  levels  of  detail 
regarding  the  atomic,  chemical,  structural,  and  mechanical  interac¬ 
tions  within  the  cell  as  are  related  to  its  composition  and  function. 
These  tools  can  provide  unique  opportunities  to  better  define  the 
system  and  conditions,  as  well  as  provide  an  effective  means  for 
validation. 

With  ongoing  experimental  and  modeling/simulation  advance¬ 
ments,  the  community  seems  poised  to  translate  individual  insights 
gained  from  systematic  studies  to  more  of  a  robust  approach  which 
considering  a  wider  breadth  of  scales,  processes,  and  interactions 
within  the  system.  This  type  of  capability  can  present  both  unprece¬ 
dented  challenges  and  opportunities  to  systematically  understand 
the  SOFC  and  transition  from  more  Edisonian  approaches  to  the 
methodical  engineering  of  advanced  materials,  structures,  and 
components  that  form  these  multiscaled,  multifunctional,  hetero- 
generous  functional  materials  (HeteroFoaMs). 
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